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[i] This article addresses the correction for aerosol effects in near-simultaneous multi- 
angle observations acquired by the Compact Reconnaissance Imaging Spectrometer for 
Mars (CRISM) aboard the Mars Reconnaissance Orbiter. In the targeted mode, CRISM 
senses the surface of Mars using 1 1 viewing angles, which allow it to provide unique 
information on the scattering properties of surface materials. In order to retrieve these data, 
however, appropriate strategies must be used to compensate the signal sensed by CRISM 
for aerosol contribution. This correction is particularly challenging as the photometric 
curve of these suspended particles is often correlated with the also anisotropic photometric 
curve of materials at the surface. This article puts forward an innovative radiative transfer- 
based method named Multi-angle Approach for Retrieval of Surface Reflectance from 
CRISM Observations (MARS-ReCO). The proposed method retrieves photometric curves 
of surface materials in reflectance units after removing aerosol contribution. MARS-ReCO 
represents a substantial improvement regarding previous techniques as it takes into 
consideration the anisotropy of the surface, thus providing more realistic surface products. 

Furthermore, MARS-ReCO is fast and provides error bars on the retrieved surface 
reflectance. The validity and accuracy of MARS-ReCO is explored in a sensitivity analysis 
based on realistic synthetic data. According to experiments, MARS-ReCO provides 
accurate results (up to 1 0% reflectance error) under favorable acquisition conditions. In 
the companion article, photometric properties of Martian materials are retrieved using 
MARS-ReCO and validated using in situ measurements acquired during the Mars 
Exploration Rovers mission. 
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1. Introduction 

[2] Imaging spectroscopy is a key remote sensing technique 
to study the composition, the mineralogy, and the physical 
state of planetary surfaces. A new generation of imaging spec- 
trometers has recently emerged in the field of space explora- 
tion by incorporating an angular dimension of measurement. 
Multi-angle imaging spectroscopy is conceived to provide a 
more accurate characterization of planetary materials from 
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orbit and a higher success in separating the signals coming 
from the atmosphere and the surface [. Martonchik et al., 
2009], The Compact Reconnaissance Imaging Spectrometer 
for Mars (CRISM) aboard the Mars Reconnaissance Orbiter 
(MRO) is a visible/short-wave infrared hyperspectral camera 
that can operate systematically in multi-angle mode from 
space [Murchie et al., 2007]. In the targeted mode, CRISM 
measures radiative energy coming from the planet Mars within 
0.36-3.92 microns using 544 spectral bands. As shown in 
Figure 1, CRISM targeted observations are composed of a 
quasi-nadir hyperspectral image at a spatial resolution of up 
to 1 8 m/pixel accompanied by ten 1 0 x spatially binned hyper- 
spectral images which are captured before and after the central 
scan. This sequence of off-nadir scans called emission phase 
function (EPF) encompasses view zenith angles up to 70° 
and two modes of relative azimuth, corresponding to the in- 
bound and outbound portions of the satellite flyby. CRISM 
targeted observations can sense a given Martian site using dif- 
ferent atmospheric paths and thus have proved to be of great 
value in atmospheric studies [Smith et al., 2009; Wolff et al., 
2009]. In addition, this multi-angle data may be exploited to 
separate atmospheric and surface contributions when analyzed 
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Figure 1 . Schema of a CRISM targeted observation. Only 
three images out of 1 1 are shown. The central scan at full spa- 
tial resolution is drawn in green, while the most extreme scans 
of the EPF are pictured in blue. Dash-dotted and dashed lines 
depict the first and last scanned lines of each single image, 
respectively. All scans are acquired as a combination of the 
gimbal and the motion of MRO which results in a VZA range 
of [0°,30°] for the quasi-nadir central scan. The two 
modes in relative azimuth (</Ji n +<p 0 ut~ 180°) are shown. 


using a radiative transfer model accounting for surface scatter- 
ing and atmospheric extinction [. Murchie et al., 2007]. The 
combination of an adequate atmospheric correction and the 
multi-angle capabilities of CRISM may make possible the 
retrieval of photometric properties of the surface according 
to observation geometry and wavelength. These surface pro- 
ducts are essential to characterize the physical state of Martian 
materials as well as to distinguish between different types of 
terrains [. Johnson et al., 2006a, 2006b; Jehl et al., 2008]. 

[3] The atmosphere of Mars is quite thin and mainly com- 
posed of carbon dioxide (C0 2 ) gas and suspended mineral 
and water ice particles named aerosols [Barlow, 2008]. 
These components represent an obstacle to remotely sensing 
the Martian surface as extinction of solar irradiance in the 
visible/short-wave infrared range happens by gaseous ab- 
sorption and aerosol scattering and absorption. The model- 
ing of the radiation reaching the top-of-atmosphere (TO A) 
level therefore demands an adequate characterization of both 
components [ Clancy and Lee, 1991; Drossart et al., 1991; 
Erard et ah, 1994; Ockert-Bell et al., 1997; Tomasko 
et al., 1999; Korablev et al., 2005; Wolff et al., 2009]. Atmo- 
spheric correction algorithms compensate the sensed TO A 
radiances for atmospheric contribution through the inversion 
of such atmosphere/surface models in order to retrieve sur- 
face quantities. One of the main hurdles to achieve an accu- 
rate correction is that aerosols and materials at the surface 
scatter light anisotropically depending on illumination and 
viewing conditions and are rarely spectrally distinct. Re- 
trieval of accurate estimates of surface reflectance therefore 
requires the consideration, in the radiative transfer model, 
of the scattering phase function of aerosols as well as the 
bidirectional reflectance distribution function (BRDF) of 
the surface. In this context, the availability of multiple 
TOA measurements with different angular configurations, 
together with the implementation of appropriate techniques, 
becomes necessary to solve as accurately as possible the 
inverse problem [Martonchik, 1994], In practice, atmospheric 
correction strategies are devised to retrieve the surface BRDF 
assuming an aerosol phase function. 


[4] While correction for atmospheric gases is generally 
possible for non-icy surfaces [Langevin et al., 2005; 
McGuire et al., 2009], compensation for aerosol contribu- 
tion represents a major challenge. Aerosol content on Mars 
is very variable as these particles are highly mobilized dur- 
ing storms which range from local events to global storms 
covering the entire planet. The variability of aerosol content 
in time and space is generally overcome by estimating, for 
each spacecraft orbital path, the local atmospheric opacity 
or the aerosol optical depth (AOD). For instance, Vincendon 
et al. [2007] propose to retrieve surface reflectance from 
observations taken by the hyperspectral imager “Observa- 
toire pour la Mineralogie, TEau, les Glaces, et TActivite” 
(OMEGA) aboard the Mars Express orbiter after quantifying 
the contribution of atmospheric aerosols using a Monte 
Carlo method. In that study, the need of multi-angle data is 
satisfied by processing several nadir observations acquired 
over the same target at different times, resulting in different 
solar zenith angles, and assuming that the AOD remains 
constant along the time span (up to several days). Besides 
this supposition, the surface is assumed to be characterized 
by a wavelength-dependent angular-independent single re- 
flectance value called Lambertian albedo. The Lambertian 
assumption supposes that variations of TOA radiance with 
geometry are exclusively related to aerosols. This hypothesis 
is also adopted by McGuire et al. [2008] to process the cen- 
tral scan of CRISM targeted observations by a radiative 
transfer approach that computes surface spectra in Lamber- 
tian albedo units after correction for gaseous and aerosol 
contributions. Alternatively, Brown and Wolff [2009] pro- 
pose to exploit the multi-angle information enclosed in 
CRISM targeted observations by modeling the remotely 
sensed signal at a single wavelength. Three parameters 
(surface Lambertian albedo, mineral aerosol opacity, and 
water ice aerosol opacity) are iteratively adjusted to fit the 
multi-angle CRISM data. Similarly, a Lambertian surface 
is also considered by the radiative transfer-based procedure 
developed by Wiseman et al. [2012] to retrieve atmospheri- 
cally corrected CRISM surface spectra. 

[5] In planetary remote sensing, the Lambertian assump- 
tion is generally adopted by atmospheric correction techni- 
ques, even when multi-angle measurements are available. 
This supposition substantially simplifies the radiative trans- 
fer modeling and reduces the size of the look-up table used 
in the inversion. The Lambertian hypothesis seems intui- 
tively appropriate in some cases (e.g., turbid atmosphere). 
However, it has been proved that non-Lambertian, or aniso- 
tropic, scattering properties play a significant role for most 
mineral and icy surfaces [de Grenier and Pinet, 1995; 
Pinet and Rosemberg, 2001; Johnson et al., 2006a, 2006b; 
Lyapustin et al., 2010]. As a consequence, the adoption of a 
Lambertian assumption can create significant angle-dependent 
biases in derived surface reflectances [Lyapustin, 1999]. To 
our knowledge, Cull et al. [2010] are pioneers in performing 
surface retrievals from CRISM targeted observations consid- 
ering a non-Lambertian surface. In that work, the inversion 
for the surface properties is achieved using a Hapke’s model 
for the surface BRDF and a radiative transfer-simulated 
look-up table that stores TOA spectra corresponding to multi- 
ple atmospheric situations. Nevertheless, the algorithm of Cull 
et al. [2010] is restricted to a few CRISM targeted observa- 
tions as it assumes that the sensed surface has similar 
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photometric properties than the materials characterized by the 
Mars Exploration Rovers (MER) mission at the Gusev crater 
of Mars. 

[6] In the field of Earth remote sensing, some efforts have 
been done towards the mitigation or the elimination of the 
Lambertian assumption. Guanter et al. [2005] address the 
atmospheric correction of multi-angle observations acquired 
by the Compact High Resolution Imaging Spectrometer 
(CHRIS) aboard the PRoject for On-Board Autonomy 
(PROBA) under a Lambertian constraint with relaxation. 
An iterative strategy for surface retrieval is defined by 
assuming a Lambertian surface on the first iteration and 
injecting the inferred BRDF into the radiative transfer-based 
inversion on subsequent iterations until convergence is 
reached. Although this method improves the quality of the 
finally estimated BRDF, it may become problematic for 
highly anisotropic surfaces for which the initial assumption 
may greatly impact the succeeding retrievals. The incorpora- 
tion of the BRDF into the atmospheric correction of 
remotely sensed data is also considered for processing 
multi-temporal data collected by the MODerate resolution 
Imaging Spectroradiometer (MODIS) aboard the Terra and 
Aqua spacecrafts [Schaaf et al., 2002] and the Multi-angle 
Imaging SpectroRadiometer (MISR) aboard the Terra space- 
craft [Diner et al., 2005]. Recently, surface retrievals are car- 
ried out for series of multi-temporal images acquired by 
MODIS by the Multi-angle Implementation of Atmospheric 
Correction (MAIAC) algorithm Lyapustin et al. [2011a, 
2011b, 2012], Contrary to the previous methods, MAIAC 
solves the aerosol/surface coupling problem without strong 
reductionist assumptions using an innovative radiative 
transfer-based formulation of the TOA radiance. The mathe- 
matical properties of this formulation combined with the con- 
sideration of an advantageous surface BRDF model enable a 
fast and efficient inversion for the surface reflectance. 

[ 7 ] To our knowledge, a method to compensate CRISM 
targeted observations for aerosol contribution considering 
an unknown non-Lambertian surface is not available in the 
literature. This article puts forward such a method through 
an innovative radiative transfer-based strategy that inherits 
some basis from the MAIAC algorithm. The method under 
the name of Multi-angle Approach for Retrieval of Surface 
Reflectance from CRISM Observations (hereafter referred 
to as MARS-ReCO) is proposed to retrieve surface reflec- 
tance from CRISM targeted observations. The major innova- 
tion is the consideration of an anisotropic surface, thanks to 
an accurate quasi-linear parametrization of the TOA radi- 
ance. In addition, the algorithm is endowed with a formalism 
integrating several sources of uncertainty into the inversion 
process and propagating them to the solution. MARS-ReCO 
transforms a CRISM targeted observation, originally in units 
of TOA radiance, into a set of photometric curves with error 
bars in units of surface reflectance depending on acquisition 
geometry. MARS-ReCO addresses the correction for mineral 
aerosols exclusively, thus not dealing with water ice aerosols 
nor atmospheric gases. The mineral aerosol particle size 
distribution and refractive index, as well as the AOD of each 
observation, must be provided beforehand. 

[8] The article is organized as follows. The MARS-ReCO 
approach is detailed in Section 2 as well as the preprocessing 
of CRISM targeted observations. Section 3 describes the tests 
conducted to evaluate the performance of MARS-ReCO on 


CRISM-like synthetic data. Eventually, Section 4 discusses 
the benefits and limitations of the proposed approach. The ca- 
pabilities of MARS-ReCO are tested on real CRISM data and 
validated against MER in situ measurements in the companion 
article [Fernando et al., 2013]. Note that the present work as 
well as the companion article restricts the use of MARS-ReCO 
to spectral bands without gaseous absorption. 

2, Methods 

[ 9 ] MARS-ReCO compensates CRISM targeted observa- 
tions for mineral aerosol effects. The proposed method inherits 
the basis from the method named MAIAC [Lyapustin et al, 
2012] while adding some functionalities. Contrary to 
MAIAC, which works with multi-temporal series of images, 
MARS-ReCO is devised to process near-simultaneous multi- 
angle CRISM observations. Furthermore, MARS-ReCO takes 
care of propagating several sources of errors to the end 
products. Before applying MARS-ReCO, CRISM targeted 
observations are transformed into appropriate products for 
multi-angular data processing. 

2.1. Preliminaries 

2.1.1. Integrated Multi-angle Product 

[ 10 ] The 1 1 hyperspectral images forming a CRISM tar- 
geted observation are released individually in units of HF, 
the ratio of measured intensity to solar flux. Note that the de- 
pendency on wavelength of all radiometric quantities in this 
article is omitted for simplicity. In order to facilitate the 
simultaneous processing of the multi-angle information 
associated to each terrain unit, the spectra corresponding to 
the 10 scans are rearranged in a common geographical grid 
of super-pixels (i.e., terrain units sensed at more than one 
geometry) to create a single data set named SPC cube (for 
SpectroPhotometric Curve). This is done after degrading 
the spatial resolution of the central scan in coherence with 
the properties of the EPF. SPC cubes are multi-angle inte- 
grated products which facilitate the access to photometric 
curves. A photometric curve is a series of radiometric mea- 
surements at a given wavelength acquired from a specific 
super-pixel with different angular configurations. The stack- 
ing of several photometric curves corresponding to consecu- 
tive wavelengths forms a spectrophotometric set. Figure 2 
(left) illustrates the typical degree of overlapping by projecting 
the footprints of the 1 1 hyperspectral images onto the common 
geographical space. Note that the overlap is non-optimal, as 
typically only ~30% of all super-pixels are sensed by four or 
more angular measurements. Figure 2(center) shows the regu- 
lar grid that is defined over the resulting mosaic such that each 
element (i.e., a super-pixel) corresponds to a single terrain unit 
typically measuring -450 x450 m 2 . This value takes into 
account the decrease in spatial resolution at the edges of the 
central and EPF images. Under this configuration, each band 
of a SPC cube contains about 2300 photometric curves, many 
of them made of less than four measurements. Note that a 
higher spatial resolution can be achieved only by tossing out 
the largest EPF super-pixels. Figure 2(right) details the struc- 
ture of a SPC cube in which each spectrum of a targeted 
observation is arranged according to its angular configuration 
(up to 1 1) along the x axis and according to its spatial location, 
or super-pixel, along the y axis. For simplicity, super-pixels 
are arranged according to the number of available 
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Figure 2. (left) Footprints of the 1 1 scans of a CRISM observation in a common geographical space. 
Note the different extent and shape of each image. In this example, the spatial resolution of the central scan 
has not been degraded for clarity, (center) Schema of the approximative overlaying grid of super-pixels. 
Each cross indicates the center of the associated super-pixel, while its color is coded according to the num- 
ber of available angular measurements, (right) Schema of a SPC cube. Each row stores the information 
belonging to a same super-pixel. Photometric curves follow the same color code that in Figure 2 (center). 
Blank positions mean absence of angular measurement. 


measurements. Under this configuration, each row at a specific 
spectral band contains a photometric curve while a x-z plane 
includes a spectrophotometric set, both products being in I/F 
units. Note that MARS-ReCO processes one spectral band at 
a time. 

2.1.2. Retrieval of Aerosol Content 

[n] MARS-ReCO must be fed with an estimate of the 
AOD associated to the CRISM observation to be processed. 
In the literature, many authors address this issue with a direct 
application to CRISM data [ Lemmon et al., 2004; Wolff 
et al., 2009; Doute and Ceamanos, 2010]. Note that the 
performance of MARS-ReCO is sensitive to the accuracy of 
the AOD estimate. Furthermore, methods avoiding or mini- 
mizing the Lambertian surface assumption should be priori- 
tized for consistency with MARS-ReCO. 

2.2. Description of the MARS-ReCO Approach 

[ 12 ] MARS-ReCO is conceived to infer photometric 
curves of the Martian surface from CRISM multi-angle 
orbital imagery. The contribution of mineral aerosols to the 
remotely sensed signal is compensated based on a radiative 
transfer-based algorithm which is fed by the scattering prop- 
erties of these atmospheric particles. After fitting a model of 
the TO A signal to the CRISM measurements, the surface 
bidirectional reflectance associated to each super-pixel in a 
SPC cube is retrieved at the acquisition geometries. The 
theoretical background and the technical aspects of the pro- 
posed approach are detailed in the following. The notation 
used in this article is listed at the end of this document. 

[13] The MARS-ReCO algorithm is based on a coupled 
surface-atmosphere radiative transfer formulation of the 
TOA signal originally proposed by Lyapustin and Knyazikhin 
[2001] and adapted by Lyapustin et al. [2011a] for the 
MAIAC algorithm. The major keystone of this accurate 
semi-analytical formulation is the consideration in the radia- 
tive transfer equation of a Green’s function to model the intrin- 
sic diffuse scattering responses of the atmosphere and a 
kernel-based scattering model for the surface. 


2.2.1. Green’s Function of the Atmosphere 

[14] Many codes solving the radiative transfer of solar 
light in the surface-atmosphere interface are not suitable 
for simulating large collections of spectra because of their 
large computational requirements. The condition for the 
lower limit of the radiative transfer (i.e., the surface reflec- 
tance) is usually part of the calculation which must be per- 
formed for each combination of physical parameters related 
to the surface and the atmosphere. The Green’s function 
method allows us to circumvent these drawbacks by a partial 
decoupling of the surface and the atmosphere. Lyapustin and 
Knyazikhin [2001] proved that this method makes possible 
the analytical combination of the atmospheric reflectivity 
and transmissivity with the surface reflectance to calculate 
TOA radiances for an arbitrary aerosol optical depth and 
acquisition geometry. This is achieved at the price of a small 
degradation of accuracy coming from two simplifications. 
Details on the use of a Green’s function in the radiative 
transfer equation and the related assumptions can be found 
in Appendix A. 

2.2.2. Expression for the TOA Reflectance 

[15] The radiance reaching the CRISM instrument (L) can be 
decomposed as a sum of the atmospheric path radiance ( D ), 
and the radiance reflected by the surface before being directly 
( L s ) and diffusely (Lf) transmitted through the atmosphere, 

L(s 0 ,s) =D(s 0 ,s) +4 s (s 0 ,s)e~ Io/l/i| +Tf(y 0 ,s), (1) 

where 5 0 and s are the incidence and view directions. The 
first term of the surface-reflected radiance can be written as 

L s (so,s) “ Sn 0 e- r °^°{p(s 0 , s) + ucopfp)p 2 (p 0 )} (2) 



where nS is the extraterrestrial solar spectral irradiance, p is 
the bidirectional reflectance factor (BRF) ( p = nf where /is 
the BRDF as defined by Nicodemus et al. [1977]), Cq is the 
spherical albedo of the atmosphere, D s is the path radiance 
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incident on the surface, and p\ and p 2 are BRF-derived func- 
tions described in Appendix A. Parameter a is a multiple 
reflection factor that depends on the surface albedo q(po) 
such that a = (l — q(go)co) ■ 

[i6] The diffusely transmitted surface-reflected radiance at 
the TOA can be calculated from L s with the help of the one- 
dimensional diffuse Green’s function of the atmosphere ( G d ) 
as it is demonstrated in Appendix A, 

L d s (s 0 ,s)= [ G d (si,s)L s (s 0 ,,si)ds l . ( 3 ) 

Jn- 


invariant kernels scaled by their corresponding kernels 
weights k= {tf,k G ,k v }. While they are not perfectly orthogo- 
nal functions, they are sufficiently independent to allow stable 
recovery of the parameters for many angular sampling distri- 
butions [Lucht et al., 2000], Furthermore, the RTLS model 
has proved to be accurate in recreating many types of natural 
surface ranging from highly anisotropic snow-covered terrains 
[Lyapustin et al., 2010] to backscattering surfaces such as 
those found in vegetation [ Schaaf et al., 2002], In particular, 
the geometric kernel models the backscattering features at 
the smallest CRISM phase angles (g w 30°). 


[17] The quantity nG d is also called bidirectional upward 
diffuse transmittance of the atmosphere. The surface albedo 
is defined as a ratio of reflected and incident radiative fluxes 
at the surface. 


q(p 0 )=F u P(p 0 )/F D ™"(p 0 ), 

F Down (p 0) = F s (p 0) + F d (p 0 ) = nSp 0 e~^ + j D s (s 0 a) lids' , 
F Up Oo) = nSp 0 e J ^q 2 (ji 0 ) + J pqi[p ')D s (so,s')ds , 

92(^0)=-/ p(s 0 ,s)pds. 

71 Jn- 


(4) 

[is] These formulas give an explicit expression for the 
TOA radiance as a function of the surface BRF. Their accu- 
racy is high, usually within a few tenths of a percent 
[ Lyapustin and Knyazikhin, 2001]. In the following, we 
work with units of top-of-atmosphere reflectance which is 
defined as R=L/(Sp 0 ) = (I/F)/g 0 . Note that MARS-ReCO 
processes SPC cubes which have been previously trans- 
formed into units of top-of-atmosphere reflectance. For 
simplicity, we use the symbol R to refer to units of TOA 
reflectance in a gas-free atmosphere populated by aerosols. 


2.2.3. Surface Scattering Model 

[19] The anisotropy of the surface is taken into account 
through its BRF modeled using a semi-empirical kernel- 
based Ross-Thick Li-Sparse (RTLS) model [Lucht et al., 
2000]. The RTLS model expands the bidirectional reflec- 
tance distribution of the surface into a linear sum of terms 
(the so-called kernels) characterizing different scattering 
modes. The kernels in the RTLS model are derived from phys- 
ical theory through simplifying assumptions and approxima- 
tions. In particular, the surface bidirectional reflectance is 
decomposed as the sum of (i) a Lambertian — or isotropic — 
contribution, (ii) a geometric component modeling the diffuse 
reflection taking into account the 3D geometrical structure of 
opaque reflectors on the surface and shadowing phenomena, 
and (iii) a volume scattering contribution simulated by a col- 
lection of dispersed facets. In this way the BRF is written as 

P(ftnP,<P) =k I + k a f G (p 0 ,p,ip) + k v f v (p 0 ,p,ip), ( 5 ) 

where subscripts refer to Lambertian (L), geometric (G), and 
volumetric (V) components, and f c , f v are predefined 
geometric kernels. More details on the RTLS kernels can 
be found in Appendix B. 

[20] The RTLS model is very advantageous to shape effi- 
cient inversion algorithms based on radiative transfer [Schaaf 
et al., 2002; Lyapustin et al., 2012] as the surface bidirectional 
reflectance is characterized by a linear combination of three 


2.2.4. Expression for the TOA Reflectance Using 
the RTLS Surface Model 

[ 21 ] The substitution of equation (5) into equations (1)— (4) 
(after normalization to reflectance units and separation of the 
kernel weights) provides the following quasi-linear expres- 
sion for the TOA reflectance: 

Rfo, F <f) = R d {Po,F v) + k L F L (p a , p) 

+k G F G {p D , p, ip) + k v F v (p a ,p, <p) + R"'(p 0 ,p), 

( 6 ) 

where R° stands for the atmospheric path reflectance coming 
from photons scattered by aerosols without interacting with 
the surface. Quantities {F L ,F G f’ v }, on the one hand, and R" 1 , 
on the other hand, are multiplicative factors for the kernel 
weights k and a nonlinear term, respectively, 

F L (p 0 ,p) = {e-^+ a p 0 l F d (p 0 )/(nS)) + G^)), 

F k {p a ,p, <p) = Q ,p,<p) + ap 0 4 Di(p 0 ,/i,<p)}e- To/ W 

(o) 

+e -to Ao G\ (p 0 ,p,ip) + xp 0 'H l k (p 0 ,p,<p), 
R n \p„,p) = acopMe-^ie-^MpM+kt-G^p) (9) 
+k a G l G l (p)+k v G]}(p)}, 

where the subscript k refers to either geometric (G) or volu- 
metric ( V) kernels. Quantities D\ D},, G av , Gq, G l v , G l J , 
G\} , Hq, and H\ represent different integrals of the incident 
path radiance (D s ) and/or the atmospheric Green’s function 
0 G d ) with the RTLS kernels (fcfv )■ The integral expression 
for these functions can be found in Lyapustin and Wang 
[2005]. Contrary to R" 1 , which strongly depends on the surface 
by means of fimetions p\ and p 2 , quantities { de- 
pend only weakly through the multiple reflection factor (a). 
Equation 6 describes the TOA reflectance as an explicit 
quasi-linear function of the RTLS kernel weights, thus pro- 
vides the means for an efficient inversion. 

2.2.5. Look-up Table 

[ 22 ] As it will be explained, the retrieval of the surface ker- 
nel weights k is possible through the inversion of equation (6) 
subject to the availability of quantities R°, {F 1 ^ ,F V }, and 
R nl . A look-up table is generated to store R° and all surface- 
independent radiometric quantities allowing us to calculate 
{F l ,F <:: ,F v } and R nl . In particular, the look-up table includes 
the basic quantities R°, D s , c 0 , and G d as well as the interme- 
diate quantities D\ D\., G av , G 1 g , G l v , G g , G]},H g , and Hy. 

[23] Basic quantities are computed using the radiative 
transfer program named Discrete Ordinates Radiative Trans- 
fer Program for a Multi-layered Plane-Parallel Medium 
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(DISORT) [ Stamnes et al., 1988]. For a homogeneous aerosol 
layer with opacity t 0 , the atmospheric reflected path reflec- 
tance ( R d ) is computed considering a dark surface (surface 
albedo set to 0) and simulating the radiation at the sensor level. 
The path radiance incident on the surface (D s ) is computed 
similarly but considering this time the downward radiance at 
the lower interface. The spherical albedo of the atmosphere 
(c 0 ) is obtained after integration of its directional-hemispheri- 
cal version which is directly provided by DISORT. Eventu- 
ally, the diffuse Green’s function of the atmosphere (G d ) is 
calculated similarly to D s but reversing the direction of light 
propagation. In other words, the atmospheric layers must be 
set in reverse order and the result must be normalized by n 
[Lyapustin and Knyazikhin, 2001]. In the case of an homoge- 
neous atmosphere, the problem for the Green’s function 
becomes identical to the problem forZ^ provided the substitu- 
tion s 0 — > s. A completely dark surface is also considered to 
avoid multiple reflections between the surface and the atmo- 
sphere. The calculation of the intermediate quantities from 
the basic functions and the RTLS kernels is explained in 
Lyapustin and Wang [2005], 

[24] A homogeneous atmosphere made of a single layer of 
mineral aerosols is used to compute the basic quantities 
look-up table. Besides its simplicity, this atmospheric model 
is appropriate to compute a “universal” look-up table repre- 
senting the average Martian atmosphere. This model supposes 
that uncertainties related to mineral aerosols (i.e., AOD esti- 
mate, particle size, and refractive indexes) are negligible and 
that water clouds are not present. If the latter particles become 
significant, a stratified atmosphere should be considered. 
However, the correction for atmospheric water ice is out of 
the scope of this work. The radiative properties of the mineral 
aerosol particles (phase function and single-scattering albedo) 
are taken from the work of Wolff et al. [2009] in which 
Martian aerosols are modeled as cylindrical particles with an 
effective radius of 1.5 /im. 

[25] All radiometric quantities stored in the look-up table 
are computed for several values of t 0 and {p 0 ,g,tp} to 
encompass different scenarios regarding atmospheric condi- 
tions and acquisition geometry, respectively. The angular 
grid density of the look-up table is selected empirically 
based on a trade-off between inversion accuracy and 
required memory. The angular range encompassed by 
CRISM targeted observations is taken into account by com- 
puting the look-up table at 8 0 e [14°,81°], 6e [0°,70°], and 

[0°, 180°], with Ap 0 = A p = 0.02 and Ap = 3°. A nearest 
neighbor technique is selected to fit CRISM observations 
to the pre-computed look-up table. In this way, we avoid a 
costly interpolation in the angular triplet {/.i 0 ,p.,<p} that 
should be done for each CRISM angular measurement other- 
wise. On the other hand, the look-up table is computed for a 
set of 12 atmospheric opacities values, namely, r 0 (l pm) = 
{0,0.05,0.1,0.2,0.33,0.5,0.75,1,1.4,2.0,2.8,4.0}. The look- 
up table is not homogeneously sampled in terms of AOD 
since the contribution of a changing aerosol content to the 
remotely sensed signal is more variable for a low AOD. A 
linear interpolation is used to obtain the look-up table values 
for a specific AOD. Under this configuration, the size of the 
look-up table is low (~50 megabytes per spectral band) due 
to its independence on the surface kernel weights k. Note 
that the look-up table is computed only once and can be then 
used for any CRISM observation. 


2.2.6. Inversion Strategy for Surface Reflectance 

[26] MARS-ReCO retrieves the photometric curve in BRF 
units corresponding to each super-pixel of a SPC cube by 
inverting the TOA reflectance model in equation (6). This 
is done based on the look-up table described above. An iter- 
ative inversion strategy is proposed based on a formalism 
borrowed from Tamntola [2005], which integrates several 
sources of uncertainty in the inversion process and propa- 
gates them to the solution. 

2.2.6. 1. Basis and Objective of the Inversion 

[27] Let R c = , . . . , j be the photometric curve in 

units of TOA reflectance measured by CRISM at a given 
gas-free spectral band. The term Ng is the number of avail- 
able angular measurements, where Ag<ll in the case of 
CRISM targeted observations. Let R' be associated to a 
given super-pixel of the SPC cube which corresponds to a 
portion of Martian surface characterized by the state vector 
k= {k L ,k G ,k l }, whose elements are unknown at first. Let k JO / 
be the RTLS kernel weights that provide the best fit between 
the observed data R c and the predicted photometric curve 
R= {R u . . -,RN g }, which is built using equation (6). The 
comparison is done by means of the root mean square error 
(RMSE) as follows: 


RMSE = 


1 Ng 


N 


K 


gj=i 



( 10 ) 


[28] The solution k so/ is therefore the set of kernel weights 
that minimizes the RMSE such that k so / = arg min RMSE. 

k 

MARS-ReCO addresses the obtention of the triplet k y „/ for 
every photometric curve in a CRISM targeted observation. 
This triplet allows us to compute the surface BRF at the 
N g sampled geometries by means of the RTLS model 
(equation (5)). Also, the existing uncertainties are taken into 
account and propagated to put error bars to the retrieved 
solution k Jo/ and the associated surface reflectance values. 


2.2. 6.2. Composition of the Inversion Matrix and 
Assumptions on the First Iteration 

[29] The model relating the state vector k to the TOA reflec- 
tances R (i.e., equation (6)) is initially dependent on the sur- 
face properties due to the nonlinear term R" 1 and, to a lesser ex- 
tent, the quantities { F L f G ,F v }. As a solution, we follow 
Lyapustin et al. [2012] who propose to retrieve k yo/ through 

an iterative inversion algorithm. Let r^ = | ;-| 0 ^ , . . . j be 
a set of reduced measurements on the first iteration n = 0 where 
rj 01 = Rj — R® — Rj lt<y> . Using equation (6) and a matrix form, 
we read 




' F m 

F <?(0) 

F n°)- 

F (0) k, 

F<°) = 

*7^(0) 

j 

pO(0) 

j 

F n o) 

j 



o' 

pG(0) 

^Ng 

F no) 

r Ng 


(ID 


[ 30 ] Lyapustin et al. [2012] remove the dependence on the 
surface by assuming a black surface on the first iteration, that 
is, k (0) = {0,0,0} . This supposition, however, results in longer 
convergence times for bright surfaces such as those found on 
the high latitudes of Mars. As a consequence, MARS-ReCO 
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assumes an isotropic surface on the first iteration k ,(l) = 
, 0} , where the subscript bck corresponds to the 
CRISM geometry in which aerosols are less predominant, that 
is, the backscattering direction. This assumption is generally 
possible because the majority of CRISM targeted observations 
are split into two modes of relative azimuth (review Figure 1), 
one of them corresponding to low ip values. 

[ 31 ] This assumption allows the iterative process to start 
by computing F (0) through the evaluation of the look-up table 
at the angular triplet {fi 0j ,pj,(pj} of each measurement^ and at 
the atmospheric input t 0 . We remember that the look-up table 
is linearly interpolated according to t 0 and a nearest neighbor 
method is used in the geometric dimension. Similarly, the 
quantity Rf (t 0 ) + Rf (k (0) , t 0 ) is computed knowing t 0 
and k <0 '. 

2.2.6.3. Characterizing the Uncertainties 

[32] In order to propagate the existing uncertainties, the a 
posteriori covariance matrix of the state vector (C/ c ) is calcu- 
lated. The state vector k is considered a random variable. For 
this purpose, we use the formalism of Tarantola [2005] in 
the framework of (i) a linear model that relates k to the set 
of reduced measurements r <0) and (ii) a set of Gaussian prob- 
ability distribution functions (PDFs) pertaining to the input 
and output parameters of the inversion problem. Each PDF 
expresses different information at different steps of the 
process (i.e., evaluation of the most probable AOD and its 
related uncertainty, a priori knowledge on the state of the 
system, CRISM measurements, and the top-of-atmosphere 
reflectance model). We note two types of uncertainties on 
vector r (0) : 

1. The error on the CRISM measurements R ‘ ' , which is 
assumed to be independent on the state of the system and 
the other geometries, with a diagonal covariance matrix 
C c with elements a], ... .a 2 N . 

2. The error on the AOD estimation t 0 , which induces an 
error on R ? + R l J llt ' 0) . In this case, we also assume the 
independence regarding the state of the system and the 
other geometries. We evaluate the elements of the covari- 
ance matrix experimentally by generating a random 
series of values according to a Gaussian PDF with mean 
t 0 and variance <r Ti .. We then calculate, based on the TO A 
reflectance model in equation (6), the population of sam- 
ples RR{ to) + Rj'ik^. t 0 ) for each geometry j £ [l,Ay. 
Likewise, the covariance matrix is computed using 
a classical estimator. 

[33] The total covariance matrix for the reduced measure- 
ments is C[ 0) = C c + C^. Last but not least, we have some 
a priori information on the solution in the form of a PDF 
characterized by its mean k (0) and a covariance matrix C k 
that we suppose to be diagonal with elements of, a\, and 
Oq taking large values. 

2.2.6.4. Retrieval of Surface Reflectance and Propagation 
of Errors 

[34] In the framework of the linear model and the Gaussian 
PDFs, the most likely a posteriori state vector k (1> and the a 
posteriori covariance matrices C kp and C lp — respectively 
associated to k (1) and the model of TOA reflectance — are re- 
trieved through the explicit least-squares solution proposed 
by Tarantola [2005]: 


k<» = k<°> 

+ Q.F (0)r ^F (0) C t F (0)r + C< 0) ) ’ (r (0) - F (0) k (0) j , (12) 

C tp = Q--QF( 0 ) r (^F( 0 )CA-F^ r + C( 0) ) ‘f(°)Q., (13) 

C,-p = F< 0 >Q p F( 0 > r + C< 0) . (14) 

[35] Vector k (1) provides a refined estimate of the surface 
BRF and, consequently, of the model F (1) . These two quan- 
tities are used in the upcoming iteration n = 1. Furthermore, 
the reduced measurement vector and its associated covari- 
ance matrix are updated through the calculation of R" l( 1 * to 
give r (1) and . This iterative process is repeated for 
several iterations before stopping on the iteration m after 
MARS-ReCO decides that the triplet of weights is satisfac- 
tory (i.e., k (m) ~ k OT/ ). As it is explained below, the number 
of iterations is set automatically according to the goodness 
of the retrieved surface BRF model. 

2.2. 6.5. Refinement of Surface Solution by Iterative Method 

[36] A convergence criterion is defined to refine the surface 
solution through several iterations. However, the physical 
sense of the retrieved BRF is not always assured. Indeed, the 
retrieved BRF model can be related to an incorrect shape 
(e.g., negative BRF values at some angles) despite a high 
goodness of fit at the measurement angles. This situation is 
likely to happen due to a deficient radiometric quality, a limited 
angular sampling, or a high aerosol content. For this reason, 
the quality of the input TOA photometric curve and the output 
retrieved surface solution is assured on each iteration by a set 
of rejection criteria. These criteria (see Appendix C for details) 
prove to reject the majority of unphysical BRF solutions. 
Convergence usually occurs after four to five iterations result- 
ing in a computational time of 1 minute to process one spectral 
band of a SPC cube at 450 m/pixel on a regular computer (i.e., 
dual 2.66 GHz quad-core processor, 6 GB RAM). 

2.2. 6.6. Characterizing the Solution 

[37] After reaching convergence on iteration m, the quality 
of the solution is characterized by a series of indicators. 
First, the root mean square error is computed to express 
the adequacy between the CRISM measurements and the 
retrieved BRF model such that 


1 Ng ~ 

RMSE^ = — V (A m) - Ff {m) kR m 'l - Fj (m) k nm) - F G{m) k G{m A 


(15) 


[38] Second, we compute the a posteriori variance on the 
estimated surface reflectance for each geometry j=\,...,N g 
as the trace of the a posteriori covariance matrix tr(C pp ) 
where C flp = QC kp Q r . Tenn Q is the linear operator 

relating the state vector to the surface BRF values p = 
|pi” ! \ . . -pjvg | (computed using equation (5)) such that 


n /. G /T1 


p( m ) = Q k W j 



(16) 
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Figure 3. Schema of the sensitivity analysis based on syn- 
thetic CRISM-like data. * Aerosol properties are provided by 
Wolff et al. [2009], *Hapke’s parameters of the minerals 
found in the Gusev crater are borrowed from Johnson et al. 
[2006a], 


[39] Finally, the a posteriori standard deviation of a pre- 
dicted BRF photometric curve is 


a 


p ~ 



(17) 


[40] In the following, this parameter proves to be an accu- 
rate indicator of the quality of the retrieved surface BRF. In- 
deed, it helps to flag those surface solutions which are wrong 
but physical and thus more difficult to detect by the rejection 
and convergence criteria. 


3. Sensitivity Analysis 

[41] In this section, the capabilities of MARS-ReCO are 
tested on simulated data following the sensitivity analysis 
depicted in Figure 3. Such task is achieved at the CRISM 
spectral band centered at X = 755.3 nm where gas absorption 
is negligible. The present research also aims at studying the 
correlation of the a posteriori standard deviation a p with 
the quality of the retrieved surface reflectances. First, a 
CRISM-like synthetic data set is generated for test in 
Section 3.1 based on the optical properties of Martian aero- 
sols and the photometric properties of some surface materi- 
als. Second, Section 3.2 evaluates the performance of 
MARS-ReCO on photometric curves with different angular 
configurations, atmospheric conditions, and surface photo- 
metric properties. Third, the requirements to process CRISM 
multi-angle data with a limited angular diversity are given in 
Section 3.3. Fourth, the stability of the inversion performed 
by MARS-ReCO is studied in Section 3.4. Lastly, the 


impact on MARS-ReCO of the uncertainties of the input 
aerosol content is explored in Section 3.5. 

3.1. Synthetic Data Set 

[42] The sensitivity analysis is performed on a controlled 
synthetic data set formed by CRISM-like photometric curves 
which are simulated in I/F units using DISORT. We con- 
sider a coupled atmosphere/surface system that is fed by 
the scattering properties of Martian dust aerosols from Wolff 
et al. [2009] and by photometric properties of Martian 
minerals. The reflectances of several surface materials found 
in the Gusev crater were measured on ground by the Pancam 
instrument aboard the MER Spirit. These data were fitted by 
Johnson et al. [2006a] to a Hapke’s BRDF model with a 
two-lobed Henyey-Greenstein phase function [Hapke, 
1993], In this research, we select four mineral materials, 
which belong to two endmembers referred by Johnson 
et al. [2006a] to as (i) endmember “Soil”, which corresponds 
to typical soils being dominant at the spatial scale accessible 
to CRISM, and (ii) endmember “Red rock”, which is related 
to rocky facets with higher anisotropic photometric proper- 
ties than soils. Table 1 details two different Pancam mea- 
surements of the same endmember at 753 nm, the first one 
corresponding to a brighter sample of the endmember 
(subindex 1) and the second one belonging to a darker sam- 
ple (subindex 2). By selecting these photometrically distinct 
materials, the sensitivity analysis aims at assessing the capa- 
bilities of MARS-ReCO against the variability of Martian 
mineral surfaces. 

[43] Based on Table 1, a TOA photometric curve is simu- 
lated by DISORT at 753 nm for each combination of the 
atmospheric/angular configurations summarized in Table 2. 
Each photometric curve is made of 1 1 HF measurements at 
constant solar zenith angle and varying view zenith angle, 
thus mimicking the acquisition of targeted observations by 
CRISM. A single combination of 11 view zenith angles, 
one for each measurement, is selected for all synthetic pho- 
tometric curves as the view angle is quite constant for all 
CRISM targeted observations. By contrast, six solar zenith 
angles are considered to embrace the angular range in which 
CRISM works, going from equatorial (0 O < 60°) to polar 
observations (# 0 > 60°). Four archetypal configurations are 
selected in terms of relative azimuth according to the typical 
functioning of MRO, each one linked to a specific couple of 
ipi and </?2 values. The fourth configuration represents the ex- 
treme case when the direction of the Sun and the direction of 
the MRO flyby are orthogonal (sp\=p 2 = 90°), resulting in 
the most limited phase angle range. Note that in this extreme 
configuration half of the angular measurements are identical 
to the other half due to the symmetry in view zenith angle. 
As regards atmospheric opacity, a set of nine AOD values 


Table 1 . Photometric Parameters at 753 nm of the Four Surface Materials Considered in the Sensitivity Analysis 


Endmember \ Hapke’s Parameters at 753 nm 

CO 

e 

b 

c 

Soil-1 (Table 4c in Johnson et al. [2006a]) 

0.69 

11 

0.241 

0.478 

Soil-2 (Table 6c in Johnson et al. [2006a]) 

0.65 

12 

0.170 

0.547 

Red rock-1 (Table 4b in Johnson et al. [2006a]) 

0.83 

19 

0.450 

0.255 

Red rock-2 (Table 8b in Johnson et al. [2006a]) 

0.65 

14 

0.166 

0.801 


A Hapke’s model with a two-lobed Henyey-Greenstein phase function is used. Parameters to, 5, b, and c correspond to the single-scattering albedo, the 
macroscopic roughness, the asymmetry parameter, and the backward scattering fraction of the phase function, respectively [Hapke, 1993]. Note that a bright 
backscattering highly anisotropic surface corresponds to to — » 1, c > 0.5, and b—> 1. 
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Table 2. Angular and Atmospheric Ranges Encompassed by the 
TOA Photometric Curves of the Synthetic Data Set a 

Parameter Acquisition and Atmospheric Configurations 

6 (°) 1 config.: 70, 63.5, 57.5, 52, 46.5, 25, 46.5, 52, 57.5, 63.5, 70 

9 0 (°) 6 config.: 30; 40; 50; 60; 70; 80 

ip u <P 2 (°) 4 config.: 0, 180; 30, 150; 60, 120; 90, 90 

To 9 config.: 0; 0.1; 0.33; 0.5; 1; 1.5; 2; 2.5; 3 

“Only one configuration is considered in terms of view zenith angle as 
this angle is quite constant among the 1 1 images of all CRISM targeted 
observations. For the central scan, we choose 0= 25° as this is the typical 
average value in real CRISM targeted observations. 

at 1 fim is considered to test MARS-ReCO under clear and 
turbid conditions. The resulting synthetic data set contains 
216 configurations for each surface material, making a total 
of 864 CRISM-like photometric curves. 

[44] Lastly, noise is added to the synthetic data set to mimic 
the signal-to-noise ratio of CRISM observations. An additive 
Gaussian noise with a noise figure of 1/50 is added to each 
simulated IIF value (i.e., the standard deviation of a given 
HF measurement R ; c being equal to R ; c /50). Note that this 
noise level is higher than the one claimed by Murchie et al. 
[2007], that is, a signal-to-noise ratio of 450 dB at 750 nm. 
In this way we assess MARS-ReCO under less favorable 
and more realistic conditions. 

3.2. Study on the Acquisition Geometry, Atmospheric 
Opacity and Surface Type 

3.2.1. MARS-ReCO Performance Based on Built-in 
Indicators 

[45] Each TOA photometric curve in the synthetic data set 
is compensated for aerosol contribution by MARS-ReCO. In 
the case of convergence, the corresponding surface RTLS 
weights k vo/ are retrieved. According to the noise attributes 

of the synthetic data set, each photometric curve R c = 

. . . ,R^,. . . , R ( ^, g | is assumed to have a diagonal co- 

variance matrix C c with elements a j = (^Rj /50^ . In order 

to determinate the intrinsic limitations of MARS-ReCO, we 
first consider the ideal case in which the exact optical depth 
(the AOD values used for the data simulation) is known. 
Accordingly, parameter is set to zero. The performance 
of the surface inversion is illustrated in Figure 4 in terns of 
(i) unsuccessful retrievals (i.e., the percentage of photometric 
curves that have been discarded for inversion), (ii) RMSE, 
and (iii) standard deviation of the retrieved BRF (a p ). 

[46] Figure 4(top) illustrates the percentage of unsuccess- 
ful retrievals according to atmospheric conditions (set by 
r 0 ), angular configuration (set by 6 0 and {v?i,v? 2 })> and type 
of material. As it can be seen, the rate of unsuccessful retrie- 
vals is about 1 5% for all TOA photometric curves not sam- 
pling the solar cross-principal plane (</rt,V J 2 = {90°,90°}). 
This specific azimuthal configuration is related to the most 
limited phase angle range and does not encompass enough 
angular diversity for a satisfactory inversion. Likewise, un- 
successful retrievals increase significantly for turbid atmo- 
spheres and extreme solar zenith angle, reaching a 40% rate 
of rejection when t 0 >2 and (9 o >80°. Regarding the solar 
angle, the RTLS model becomes less accurate at extreme 


illumination conditions due to the divergence of the geomet- 
ric kernel [Lucht et al., 2000] and the higher anisotropy of 
surfaces in this angular range. The success of MARS-ReCO 
is, by contrast, less dependent on the properties of the sur- 
face materials. In particular, MARS-ReCO performs slightly 
worse when dealing with bland materials such as Soil-2, 
which may be too dark and too isotropic to be accurately 
separated from the aerosol photometric curve. 

[47] Figure 4(middle) illustrates the quality of the fit be- 
tween the synthetic data set and the TOA reflectance model 
by means of the RMSE (see equation (15)). Only the suc- 
cessfully inverted photometric curves are considered in this 
experiment. The RMSE is generally equal to a few tenths 
of 1% reflectance, and while it logically increases according 
to solar zenith angle and AOD, it decreases for reduced 
phase angle ranges. This result is, however, reasonable since 
photometric curves along the solar cross-principal plane are 
easier to fit as they do not sample the aerosol photometric 
curve spanning the solar principal plane (yq, ip 2 = {0°. 180°}). 
Contrary to the lower RMSE, the retrieved BRF model in 
this case is likely to be wrong. Therefore, the RMSE 
cannot be considered as a reliable indicator of the quality of 
MARS-ReCO as the surface/atmo sphere model can satisfacto- 
rily fit a photometric curve while providing a physically incor- 
rect solution. By contrast, the RMSE coherently increases 
according to solar zenith angle due to the combined effect of 
a stronger surface anisotropy at extreme angular configura- 
tions and the limitations of the RTLS model. Note that inaccu- 
racies at extreme angles also come from small differences 
between the Hapke’s and the RTLS models as well as 
the use of a plane-parallel radiative transfer code such as 
DISORT. Lastly, the RMSE also depends slightly on the type 
of surface, being higher for anisotropic surfaces. 

[48] Finally, Figure 4(bottom) explores the average a pos- 
teriori standard deviation of the surface photometric curves 
(a p ) as defined in equation (17). Again, only the successfully 
inverted photometric curves are considered. The standard 
deviation is generally lower than 0.1 for favorable condi- 
tions and, contrary to the previous indicators, clearly 
increases according to AOD, solar zenith angle, and limited 
phase angle range, which typically correspond to the most 
challenging scenarios. 

3.2.2. MARS-ReCO Performance Based on the Pancam 
Reference 

3.2.2. 1. BRF Error in the 11 Acquisition Geometries 

[49] The accuracy of the retrieved surface BRF is evalu- 
ated using the surface data derived from Pancam. Note that 
these data are used in the simulation of the synthetic data 
set. This experiment is done exclusively for those photomet- 
ric curves for which MARS-ReCO converges. We define the 
BRF error of a given photometric curve ( e p ) as the average 
of the difference between (i) the surface photometric curve 
retrieved by MARS-ReCO in BRF units and (ii) the same 
curve reconstructed using the Hapke’s model and Table 1 

1 oo^|pf ARS ' ReCO -pf apke | 

e P — jUp. 2-^ Hapke ' 

S ,= 1 Pj 

[50] Figure 5 shows the result of averaging the BRF error 
of all photometric curves associated to a given AOD, solar 
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Sun zenith angle Azimuthal configuration 


Figure 4. (top) Percentage of unsuccessful retrievals, (middle) RMSE, and (bottom) a posteriori standard 
deviation of the retrieved BRF for the processed synthetic data set. All parameters are averaged according 
to the AOD, the solar zenith angle, and the azimuthal configuration that correspond to each synthetic pho- 
tometric curve. The RMSE and the standard deviation are computed based only on successfully inverted 
curves. 


zenith angle, or azimuthal configuration. As it can be seen, the 
accuracy of the retrieved BRF decreases under unfavorable 
conditions, that is, high AOD, extreme solar zenith angle, 
and limited phase angle range. The average BRF error is 
lower than 20% when t 0 < 1 , 9 0 < 60° , and ip i , p 2 ^ { 90° , 90° } 
despite the unfavorable configurations of the other param- 
eters (e.g., e p > 20% when t 0 < 1 and 6 0 = 80°). By contrast, 
it is important to remark that ranges of validity with higher 
average BRF errors contain configurations with acceptable 
accuracies (e.g., for Soil-2, while the average e p is up to 
38% when t 0 = 2, the individual e p is 5.6% for the configu- 
ration when t 0 = 2 , 9 0 = 30°, and p\, p 2 = {30°, 150°}). 
Note the high correlation between the BRF error and the a 
posteriori standard deviation ( a p ) in Figure 4(bottom). 
Lastly, we remark the continuity of the error curves when 
to =1.5. Contrary to the adjacent situations (t 0 =1 and 
r 0 = 2), this particular aerosol content is not considered in 
the look-up table, and therefore, the inversion of the associ- 
ated curves is made after interpolation of the pre-computed 
atmospheric quantities (see Section 2.2.5). According to 
results, the performance of MARS-ReCO is robust against 
this interpolation. 


3.2.2.2. BRF Error in the Complete Upper Hemisphere 

[ 51 ] The capabilities of MARS-ReCO are further validated 
by evaluating the retrieved surface BRF model out of the 1 1 
acquisition geometries. This time we use a dense angular grid 
to explore the whole upper hemisphere in view zenith angle 
and relative phase angle. The evaluation is, however, restricted 
to the illumination conditions of acquisition (6 0 ) as the inves- 
tigation of other angular ranges becomes unpredictable. 
Figures 6 and 7 assess a couple of surface models that have 
been retrieved under varied acquisition conditions. The BRF 
models are evaluated at the defined angular grid and plotted 
in polar coordinates, where the radial and angular coordinates 
correspond to the view zenith angle and the relative azimuth of 
evaluation, respectively. Note that while the case illustrated 
in Figure 6 corresponds to very favorable conditions, the 
examples in Figure 7 are increasingly challenging in tenns 
of AOD, solar zenith angle, azimuthal configuration, and 
surface anisotropy. 

[52] Figure 6 illustrates the evaluation of the retrieved 
BRF model for a surface made of material Soil-1 in the ab- 
sence of atmosphere (t o = 0). The Sun position is at 6 0 = 30° 
(see yellow star) and the MRO flyby results in a relative 
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Aerosol optical thickness Sun zenith angle Azimuthal configuration 


Figure 5. BRF error in percent (e p ) between the BRF curve retrieved by MARS-ReCO and the reference 
data constructed from a Hapke's model fed by Table 1. The plotted BRF errors are the result of averaging 
each combination of (left) AOD, (center) solar zenith angle, and (right) azimuthal configuration. Only 
successful retrievals are considered. 


Reference Hapke's model 



VZA 


Retrieved RTLS model BRF Error e p 




VZA VZA 


Figure 6. Surface BRF corresponding to material Soil-1 generated using (left) the reference Flapke’s 
model and (center) the RTLS model retrieved by MARS-ReCO from a TOA photometric curve defined 
by t 0 = 0 (i.e., no atmospheric effects), 0 o = 3O°, and <P\,<P 2 = {30°, 150°}. (right) BRF error ( e p ) between 
the BRF calculated from the Flapke’s model and that from the RTLS model. Red stars show the geometry 
of the 1 1 measurements of the synthetic photometric curve. The Sun position is marked with a yellow star. 
Note that the backscattering direction is situated at the right-hand side of the plot. 


azimuthal angle configuration of < / 9 1 ,</? 2 = {30°, 150°}. Red 
stars represent the 1 1 CRISM-like measurements. This angu- 
lar configuration results in a phase angle range of acquisition 
where ge[14°,96°]. Figure 6(left) shows the surface BRF 
calculated using the Hapke’s model fed by the reference pho- 
tometric data in Table 1 and evaluated at 00 = 30°. Similarly, 
Figure 6(center) shows the surface BRF calculated using the 
RTLS model fed by the k vo/ retrieved by MARS-ReCO. 
Figure 6(right) expresses the quality of the retrieval by plotting 
the relative difference between both BRF data sets. A low 
error (e p = 0.8%) underlines the compatibility of the RTLS 
and the Hapke’s models even with a few available measure- 
ments and despite the additive noise. This result is in agree- 
ment with the low BRF standard deviation ( a p = 0.002) given 
by MARS-ReCO. 

[53] The same experiment is repeated in Figure 7a for a 
typical atmospheric opacity on Mars, that is, i o = 0.5. 
Despite the aerosol effects, the BRF error and the standard 
deviation are very low (e p =1.6% and a p = 0.004), mainly 
due to the favorable angular and atmospheric conditions. 
As it can be seen, MARS-ReCO satisfactorily retrieves a 
surface model with two scattering lobes in the backward 
and forward directions. The sole dissimilarity is observed 
in the forward direction when 6 « 80° where the differences 
between the two surface models become significant. 


[54] In Figure 7b, the aerosol contribution is severely in- 
creased (t 0 = 2) to recreate very turbid conditions. This config- 
uration results in a higher BRF error (e p = 5.6%), which is well 
correlated with a higher standard deviation (a p = 0.027). Note 
that the BRF error is relatively moderate since strong inaccu- 
racies mainly happen at high view zenith angles (errors up to 
50%), which are not sampled by the CRISM geometries. 

[ 55 ] Figure 7c repeats the same experiment with t o = 0.5 
and 0 o = 7O° (gG [28°, 130°]). Although the retrieved back- 
scattering lobe is accurate, the extreme illumination condi- 
tions result in a somewhat incorrect forward scattering 
feature. Nonetheless, similar to Figure 7b, the average BRF 
error along the 1 1 geometries is moderate (e p = 3.0%) as well 
as the associated standard deviation {a p = 0.011). Note that 
the BRF error is expected to increase for higher values of 
solar zenith angle due to limitations of both surface models. 

[ 56 ] The experiment in Figure 6 is repeated with t o = 0.1 
and pi,tp 2 = {90°, 90°}, resulting in a very limited phase 
angle range of acquisition (gG [38°, 73°]). The inversion of 
this photometric curve is one of the few that converges under 
this azimuthal configuration [Figure 4(top right)], probably 
because the rest of acquisition conditions are favorable. 
Nonetheless, Figure 7d shows that the shape of the retrieved 
surface model is strongly inaccurate (error up to 80%). This 
result comes from the restriction of the TOA measurements 
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Figure 7. From top to bottom, same as Figure 6 when (a) t o = 0.5, (b) 
(d) t 0 = 0. 1 and </?i,<p 2 = {90°, 90°}, and (e) t 0 = 0.5 and surface made 


80 53 27 0 27 53 80 

VZA 

r 0 = 2, (c) t 0 = 0.5 and 6 0 
of material Red rock-1. 


= 70°, 


to the solar cross-principal plane where the main surface 
scattering features of Soil- 1 cannot be sampled. Note, how- 
ever, that the average BRF error and the standard deviation 
are quite low (e p = 1.8% and a p = 0.003) as inaccuracies are 
found out of the 1 1 sensing geometries. This example high- 
lights the differences between the quality of the retrieved 
photometric curve, which is quite high in this case, and that 
of the complete BRF model, which is very poor. 


[57] The last experiment summarized in Figure 7e considers 
t o = 0.5 and material Red rock-1, which is related to a higher 
albedo, a higher anisotropy, and a higher surface roughness 
(Table 1). As it can be seen, the retrieved BRF model is rather 
accurate and reproduces accurately the main narrow backscat- 
tering lobe (i.e., the scattering feature in the forward direction 
of Red rock-1 is smoothed by the surface roughness 0), except 
for a slight underestimation of 7%. The smoothness of narrow 
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and strong lobes is a typical limitation of the RTLS model 
[Roujean et al., 1992]. Again, the highest BRF error is 
obtained in the forward direction for extreme viewing geome- 
tries ( 0 > 70°). This configuration is associated to a low aver- 
age BRF error and a low standard deviation (e p = 2.5% and 
<r p = 0.005), indicating the goodness of the retrieval. 

3.3. Study on Observations with Restricted 
Geometry 

[58] Up to this point MARS-ReCO has been applied to 
photometric curves made of 11 measurements. Unfortu- 
nately, this situation does not correspond to reality since less 
than 20% of the area encompassed by CRISM targeted 
observations is typically sensed by more than seven geome- 
tries (see Figure 2). The following experiments aim at asses- 
sing MARS-ReCO against a reduced number of measure- 
ments and thus a restricted geometry. 

3.3.1. Number of Angular Measurements and Phase 
Angle Range 

[59] The synthetic data set defined in Section 3.1 is itera- 
tively degraded by removing an increasing number of mea- 
surements from each photometric curve. The position of 
the removed measurements in terms of view zenith angle is 
chosen randomly. MARS-ReCO is applied to the resulting 
data set on each iteration until only three measurements are 
available. Note that the inversion for the RTLS kernel 
weights with less than three measurements becomes uncon- 
strained. The performance of MARS-ReCO is assessed 
according to the number of measurements and the phase an- 
gle range encompassed by each photometric curve. Note that 
a high number of measurements within a narrow phase angle 


range may be less appropriate than few measurements and a 
broad phase angle range [ Soitchon et al., 2011], 

[60] According to Figure 8(top left), the rate of unsuccess- 
ful retrievals remains rather acceptable except for photomet- 
ric curves with a phase angle range lower than 40° . In this 
case, photometric curves do not contain enough angular di- 
versity to separate the photometric curve of aerosols from 
that of the surface, regardless of the number of geometries. 
Note how the performance of MARS-ReCO decreases for 
a reduced number of geometries as the inversion becomes 
progressively unconstrained. However, the number of avail- 
able geometries is less crucial than the phase angle diversity 
in terms of convergence. 

[61] Figures 8(bottom left) and 8(bottom right) explore the 
a posteriori standard deviation (er p ) and the error ( e p ) of the 
retrieved BRF, respectively. Again, the phase angle range 
appears as the most critical parameter to perform an accurate 
atmospheric correction. In fact, Figure 8(bottom right) 
shows the decrease of the retrieved BRF accuracy for limited 
phase angle ranges combined with a few available geome- 
tries. These unfavorable situations can be detected when 
processing real CRISM observations based on the a poster- 
iori standard deviation. Indeed, this parameter shows a 
strong correlation with the BRF error, thus confirming its 
pertinence as an indicator of the accuracy of the retrieved 
surface. By contrast, Figure 8(top right) proves that the 
RMSE is not correlated with the accuracy of MARS-ReCO. 
For instance, contrary to the BRF error, the RMSE is mini- 
mum when only three geometries are available as photomet- 
ric curves are easier to fit in this case. Lastly, note that the 


Unsuccessful retrievals (%) 


RMSE 


x 10 



11 10 9876543 
Number of angular measurements 


11 10 9876543 
Number of angular measurements 


Figure 8. (top left) Percentage of unsuccessful retrievals, (top right) RMSE, (bottom left) a posteriori 
standard deviation of the retrieved BRF surface er p , and (bottom right) average BRF error (e p ) computed 
along the acquisition geometries according to the number of measurements and the phase angle range 
(Ag). The RMSE, BRF error, and parameter (a p ) are computed only for the successfully processed photo- 
metric curves. The y axis of the figures is defined such that the row associated to a phase angle range of 
110° encompasses the photometric curves within 100° <Ag< 120°. Results when Ag= 70° and Ag= 170° 
are not shown due to the absence of photometric curves at these configurations. 
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rather inaccurate results obtained when dealing with very 
broad phase angle ranges (larger than 140°) are misleading. 
In fact, this situation encompasses the most extreme illumi- 
nation conditions ( 9 0 > 70°). 

3.3.2. Distribution of the Phase Angle 

[ 62 ] The previous experiment does not take into account 
which phase angles are sampled by the photometric curves. 
Although broad sampled phase angle range are usually ben- 
eficial, not all ranges of a given magnitude are equally use- 
ful. Furthermore, since the end of 20 10, CRISM is acquiring 
targeted observations without inbound portion due to gimbal 
stickiness [. Murchie , 2012], These restricted range observa- 
tions sample either small phase angles (when the Sun is 
“behind” MRO) or large phase angles (when the Sun is “in 
front” of MRO) while corresponding to similar magnitudes 
of phase angle range. 

[63] In order to assess the performance of MARS-ReCO 
according to the sampled phase angles, two synthetic data 
sets are built based on the data set described in Section 3.1. 
Both data sets mimic the latest restricted range observations 
by containing photometric curves with only five measure- 
ments corresponding to a single mode of azimuth, ipi or tp 2 
(see Table 1). The first data set simulates no inbound obser- 
vations sampling small phase angles, while the second one 
favors larger phase angles. The BRF error at the acquisition 
geometries e p resulting from the inversion of these two data 
sets is explored according to the type of material in Figure 9. 
Two cases are defined to study the sampled phase angles. The 
first case computes the average BRF error of all “backward” 
photometric curves (g E [0° , 90°]), while the second case con- 
siders all “forward” curves (g E [90°, 180°]). Note that there 
are approximately 30% less photometric curves in the second 
case as the phase angle becomes quite small for ip 2 = 120° or 
90°, especially when solar zenith angle is small. 

[64] According to Figure 9, the predominance of aerosol 
contribution results in larger BRF errors for “forward” photo- 
metric curves. Also, the BRF error is larger where there is less 
signal coming from the surface as it happens for backscattering 
surface materials when g E [90°, 180°]. Contrarily, the BRF 



g<90 deg g>90 deg 

"Inbound-only" photometric curves 


Figure 9. Average BRF error at the acquisition geometries 
(e p ) for a set of photometric curves without the inbound por- 
tion found in CRISM targeted observations. Two classes are 
defined according to the distribution of the phase angle, 
those sensing the surface of Mars exclusively (i) in the 
backward direction (g E [0°,90°]) and those (ii) in the 
forward direction (g 6 [90°, 180°]). 


error related to Red rock-1 is larger when g E [0°, 90° ] as 
this material is forward scattering. The BRF error decreases 
for large reflectances and sharp scattering lobes (e.g., Red 
rock-1) as the surface signal is easier to retrieve in this case. 
In opposite, flat scattering lobes result in higher BRF errors 
(e.g., Soil-2 obtains a higher error than Soil-1 since the latter 
has a sharper lobe) even for highly anisotropic materials 
(e.g., Red rock-2). In conclusion, while the BRF error for 
no inbound observations is somewhat higher than for com- 
plete observations (see Figure 5), it is still reasonable 
(e p «20%). This result underlines that MARS-ReCO is 
appropriate to process restricted range CRISM observations 
provided that five angular measurements are available. In 
particular, backscattering surfaces are prone to be retrieved 
more accurately when the outbound images have the Sun in 
the back (i.e., small phase angles), that is, when MRO in is 
the north (south) hemisphere in its descending (ascending) 
node. Finally, it should be remarked that the retrieved BRF 
model is likely to be less accurate out of the acquisition 
geometries than for complete targeted observations. 

3.4. Study on the Stability of the Surface Solution 

[65] This experiment explores the robustness of MARS- 
ReCO by assessing the stability of the retrieved surface solu- 
tion for a set of similar photometric curves. With this aim, 
MARS-ReCO is run on a set of photometric curves deriving 
from a particular configuration defined by material Soil-1, 
t o = 0.5, 0 o = 3O°, and (pi,<P 2 = {30°, 150°}. Figure 10 shows 
the evolution of the retrieved RTLS kernel weights accord- 
ing to a varying AOD, solar zenith angle, and azimuthal 
configuration. In Figure 10(left), for instance, the AOD asso- 
ciated to the set of photometric curves varies in contrast to 
the solar zenith angle and the azimuthal configuration, which 
remain at their initial values. According to results, the surface 
BRF models provided by MARS-ReCO for material Soil-1 
are quite alike despite the variability encompassed by the 
data set of study. Only extreme configurations such as ipi,ip 2 = 
{90°, 90°} or 9 0 = 80° result in a significantly different RTLS 
triplet k v „/. As for the aerosol content, Figure 10(left) shows 
that an increasing AOD results in a slight variation of all RTLS 
kernel weights (note that the geometric kernel function F° is 
much larger than the volumetric one F 1 , and therefore, the 
associated weights k G are respectively smaller). The increase 
in k 1 is compensated by the decrease in k G , so the anisotropy 
of the surface is maintained. This result underlines that two 
slightly different combinations of the RTLS weights can pro- 
vide two satisfactory surface fits at the acquisition geometries. 
Note, however, that when t 0 > 2 the retrieved BRF models are 
certainly incorrect. 

3.5. Study on the Aerosol Content Uncertainty 

[66] Up to now MARS-ReCO has been fed with the exact 
aerosol content that is used for the data simulation. In reality, 
however, the actual optical depth of a given observation can 
be known only with a given uncertainty. Clancy et al. [2003] 
find a ±0.05 uncertainty for retrieved optical depths between 
0.20 and 0.50 using EPF observations taken by the Thennal 
Emission Spectrometer aboard Mars Global Surveyor. 
Lemmon et al. [2004] determine typical uncertainties 
between ±0.02 and ±0.04 when measuring the AOD with 
the MER rovers. Vincendon et al. [2009] find uncertainties 
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Figure 10. RTLS kernel weights retrieved from a set of photometric curves initially defined by material 
Soil-1, t o = 0.5, 0 o = 3O°, and ip\,tp 2 = {30°, 150°}. The AOD, solar zenith angle, and azimuthal configu- 
ration vary in the left, center, and right figures, respectively. 


of ±10% for moderate aerosols loading (t o (1/hm) = 0.5) and 
of ±20% for high aerosols contribution (t 0 ( 1 pm) = 1) using 
OMEGA observations. Similarly, Wolff et ai [2009] obtain 
uncertainties of ±10-20% for aerosol depths measured with 
CRISM EPF data. 

[67] The last experiment of the sensitivity analysis investi- 
gates the impact of these uncertainties on the MARS-ReCO 
capabilities. We now consider that the actual optical depth is 
affected by an unknown additive bias such that t 0 = 
t 0 ± b Zo . According to the studies cited above, the bias af- 
fecting aerosol content is modeled with a normal distribution 
that introduces an average bias of ± 1 5% (standard deviation 
of 5%). In this way, a vector of 216 different biases b Zo — 
negative and positive with different magnitudes — is gener- 
ated, one for each synthetic photometric curve. The suite of 
inversions performed in Section 3.2.1 is repeated using the 
biased optical depths (t' 0 ). The error on the latter estimates 
is considered by setting the parameter characterizing the 

AOD uncertainty to <r^ o = ^0.15 * t 0 ^) (see Section 2.2.6). 

[68] Figure 1 1 illustrates the inversion quality by means 
of the typical MARS-ReCO indicators. Only the dependence 
on aerosol optical depth is shown as the relation between 
the inversion quality, and all acquisition parameters is very 
similar to that observed for the ideal case (see Figures 4 
and 5). In fact, only the average value of the quality indica- 
tors suffers significant variations. In particular, the RMSE 


and the BRF error (e p ) undergo an increase of 30% and 
9%, respectively. This decrease in the quality of the retrieved 
surface reflectance results in an increase of the evaluated 
uncertainty of the surface BRF (cr p ) by a factor 2. Higher 
inaccuracies for r 0 = 2 are explained by the higher bias 
that affects the photometric curves typically related to 
the greatest errors (e.g., d 0 = 80°) in this AOD interval. Like- 
wise, remember that the CRISM-like data are affected by 
random noise. Lastly, the rate of unsuccessful retrievals is 
not shown for this experiment as results are very similar 
to those in Figure 4. Although the rate of successful retrie- 
vals decreases due to the error on the AOD estimates, the 
relaxation introduced in the inversion by considering ^ 0 
helps MARS-ReCO fitting the TOA photometric curves. 

4. Discussion 

[69] According to the sensitivity analysis described in the 
previous section, MARS-ReCO proves to be appropriate to 
process real CRISM targeted observations. Surface reflectance 
curves provided by MARS-ReCO are physically plausible and 
mostly accurate. The few inaccurate surface solutions can be 
detected by examining the a posteriori standard deviation of 
the estimated BRF ( o p ) provided by MARS-ReCO. This 
output has proved to be a satisfactory indicator of the quality 
of the retrieved surface reflectance according to its high corre- 
lation with the BRF error. Note that the latter indicator is 




Aerosol optical thickness 


Figure 11. RMSE, a posteriori standard deviation of the retrieved BRF, and BRF error when a 10% un- 
certainty is considered for the aerosol content estimate. Note the different vertical axis range with respect 
to Figures 4 and 5. 
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accessible only if ground truth data are available. Based on the 
sensitivity analysis, this section establishes the range of condi- 
tions over which MARS-ReCO is able to confidently retrieve 
surface properties from CRISM targeted observations. 

[ 70 ] Table 3 summarizes the main results obtained by the 
sensitivity analysis. In particular, retrievals of surface reflec- 
tance are possible (unsuccessful retrievals lower than 5%) 
and relatively accurate (BRF error lower than 20%) when 
TOA photometric curves are related to (i) an appropriate 
angular configuration set by a moderate solar zenith angle 
( 6 0 <70°) and an azimuthal configuration ensuring a broad 
phase angle range (Ag e [40°, 140°]), and (ii) a surface sig- 
nal which is significant compared to the aerosol contribu- 
tion, that is, a moderate AOD (t 0 < 1 .5). Accuracy can be 
increased (e p < 10%) by restricting the validity ranges to 
0o <60° and r o <0.9. Furthermore, results show that surfaces 
combining a low albedo and a low anisotropy (e.g., Soil-2 
with co = 0.65 and c = 0.17) lead to slightly worse retrievals 
(accuracy decrease of 10%). The previous range of condi- 
tions is subject to the availability of the aerosol content over 
CRISM observations. An average decrease of 9% of the 
quality of the retrieved surface reflectance must be expected 
for AOD estimates with a 15% error. 

[ 71 ] In practice, the range of conditions providing inaccu- 
rate surface solutions is limited. First, the inversion of CRISM 
observations with restricted phase angle range (e.g., yq, <p 2 = 
{90°, 90°}) is discarded by the criteria of MARS-ReCO 
(90% of rejection in this case). In this matter, experiments 
underline that the reflectance accuracy remains acceptable 
( e p ;S 20%) even against a reduced number of angular mea- 
surements (down to 5-6 geometries) as long as a broad phase 
angle range is available (Ag G [40°, 140°]). Configurations 
out of these bounds make the separation between aerosol 
and surface signals not possible due to the lack of angular 
diversity, when Ag < 40°, or the scarcity of available geome- 
tries, when Ag> 140°. Second, opacities on Mars are often 
rather low (r 0 < 1 long-ward of one micron [Smith, 2009]). 
However, precaution must be taken when processing CRISM 
observations from the high latitudes of Mars as they can be 
acquired with solar zenith angles around 70°. Besides the 
limitations of the RTLS surface model at these angles, the 
plane-parallel approximation used by DISORT may affect 
the surface retrievals. In this situation (and other unfavorable 
cases), the a posteriori standard deviation of the estimated 
BRF (o p ) must be inspected in order to detect incorrect 
surface retrievals. Experiments prove that surface solutions 
with Op 55 0.08 are related to errors greater than 20%. 

[ 72 ] It is important to remark that despite a low standard 
deviation (a p < 0.08), the retrieved BRF surface model must 
be handled with care when it is evaluated out of the CRISM 
acquisition geometries (Figure 7d). Also, MARS-ReCO has 


proved to be stable (variability of -10% of the retrieved sur- 
face model) when dealing with a group of adjacent CRISM- 
like super-pixels belonging to the same surface material but 
observed with slightly different geometric configurations. 

5. Conclusions 

[ 73 ] A radiative transfer-based algorithm for atmospheric 
correction is put forward to retrieve surface reflectance from 
CRISM multi-angle imagery. The method referred to as 
MARS-ReCO transforms a set of gas-free TOA radiances, 
which correspond to the same terrain unit observed at differ- 
ent geometries, into a photometric curve in units of surface 
reflectance. MARS-ReCO represents a substantial improve- 
ment regarding the state of the art in atmospheric correction 
of Martian data as it considers a non-Lambertian surface 
along with the photometric curve of the aerosols. The latter 
data are an input of the proposed technique as well as the 
aerosol optical depth. Although being based on a radiative 
transfer model, MARS-ReCO is very fast as it is based on 
an efficient formulation of the TOA signal and a linear 
model for the surface reflectivity. In addition, MARS-ReCO 
integrates a statistical formalism that propagates several 
uncertainties to the solution, thus providing an accurate indica- 
tor of the quality of the retrieved surface reflectance. These 
statements are confirmed by the sensitivity analysis that is pre- 
sented in this article and that has proven the validity, the accu- 
racy, and the stability of the surface solutions retrieved from 
realistic synthetic data (CRISM-like geometric configuration 
and noise) under realistic conditions (Martian optical depths 
associated to the typical uncertainties provided by aerosol 
retrieval methods). This is further proved in the companion 
article [Fernando et al., 2013] where MARS-ReCO processes 
real CRISM observations and retrieves surface reflectance 
curves, which are highly accurate (maximum error of 5%) 
when compared to in situ measurements. 

[ 74 ] MARS-ReCO is subject to some limitations when 
processing CRISM multi-angle data acquired under unfavor- 
able acquisition conditions, that is, turbid atmospheres, 
extreme illumination conditions, and restricted phase angle 
ranges. The uncertainties on the surface reflectance retrieved 
under these situations are, however, quantified by the a 
posteriori standard deviation provided by MARS-ReCO. It 
is important to note that the potential of MARS-ReCO to 
map the photometric properties of the surface of Mars is sub- 
ject to the degree of overlap among the individual images 
composing a CRISM targeted observation. In this matter, 
MARS-ReCO has proved to be robust against observations 
with restricted geometry such as those acquired by CRISM 
from late 2010. Furthermore, restricted geometry diversity 
can be overcome by combining different CRISM observations 


Table 3. Performance of MARS-ReCO in Terms of Unsuccessful Retrievals and BRF Error (e p ) According to Acquisition Configuration 


Convergence and BRF Error 

Optical Depth 

Solar Zenith Angle 

Phase Angle Range 

Unsuccessful retrievals < 5% 

to <2.0 

e Q < 70° 

Ag> 40° 

Unsuccessful retrievals > 5% 

t o >2.0 

$ 0 >70° 

Ag < 40° 

e p < 10%, (<j p < 0.04) 

t 0 <0.9 

#o<60° 

Ag < 100°“ 

e p < 20%, ((7 p £ 0.08) 

T 0 < 1.5 

e 0 < 70° 

Age [100°, 140°] 

e p > 20%, ( a p & 0.08) 

To > 1.5 

0 o >7O° 

Ag > 140° 


“The few photometric curves satisfying Ag < 40° that pass the inversion criteria (see in Appendix C) lead to physically plausible surface reflectance 
curves (low e p ) but strongly incorrect BRF models (see Figure 7d). 
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of the same target acquired under different angular configura- 
tions. Eventually, note that this article assesses the perfor- 
mance of MARS-ReCO on gas-free wavelengths. In order to 
process the whole CRISM spectral range, gaseous contribu- 
tion must be previously compensated [ McGuire et al., 2009; 
Doute, 2009]. In this case, testing must be carried out to eval- 
uate the impact of the uncertainties coming from a faulty gas 
correction on the retrieved surface reflectance. 

[75] One major advantage of MARS-ReCO is that its 
look-up table can be easily adapted to accommodate differ- 
ent scenarios. In this article, the look-up table considers only 
mineral aerosols, making MARS-ReCO suitable to process 
most CRISM observations acquired over the equatorial 
regions of Mars. Future versions of the look-up table may, 
however, consider other atmospheric situations such as the 
presence of water ice aerosols, mineral aerosols with differ- 
ent grain size, or even gases. Further testing should be per- 
formed to test the performance of MARS-ReCO in this case. 
Likewise, limitations found at extreme zenith angles may be 
mitigated using a look-up table calculated by a radiative 
transfer code considering the sphericity of Mars. 

[76] Further research will be conducted on the exploitation 
of the complete surface reflectance model retrieved by 
MARS-ReCO to calculate directional-hemispherical albedos 
of the surface as it is done with in situ measurements in Bell 
et al. [2008]. Also, we will address the processing of the 
whole spectral dimension of CRISM. The resulting increase 
of the input data, together with the consideration of appro- 
priate techniques, may allow us to retrieve simultaneously 
the surface BRF and the aerosol optical depth as it is done 
by Lyapustin et al. [2011b]. 


Appendix A: Green’s Function Method for the 
Radiative Transfer Problem 

[77] Lyapustin and Knyazikhin [200 1 ] address the use of the 
Green’s function method for the radiative transfer problem. 
The surface is considered to be non-Lambertian and spatially 
homogeneous, and the atmosphere is vertically stratified with 
a plan parallel geometry and t 0 being the integrated optical 
depth. The atmosphere is illuminated at the TOA by a quasi 
collimated solar beam of flux nS with incident direction 
So = (dofi). The direction of radiation propagation within the 
atmosphere is noted s = (9,p). The vertical axis z points down- 
ward so that downward directions of propagation (p > 0) are 
indicated by a plus sign (+) while the upward direction (p < 0) 
corresponds to a minus sign (— ). The wavelength of measure- 
ment of the sensor X is omitted for clarity. The method is based 
on the additive properties of the radiative transfer equation and 
on three ideas. 

1 . First, the total radiance at a given level of optical depth t 

and direction of propagation s is decomposed into the ra- 

diance corresponding to the photons that never interacted 

with the surface and a term corresponding to the radiance 
carried by all other photons 

L(z;s 0 ,s) = L>(z;s 0 ,s) + J(z;s 0 ,s), (Al) 

where the first term depends only on the properties of the 

atmosphere and obeys the radiative transfer equation with 


the boundary conditions D+( 0) = 0 and D(z 0 ) = 0. 

2. Second, the determination of the radiative response of the 
atmosphere to a given field of upward radiance at the 
surface-atmosphere interface T_(t 0 ,/) requires the calcu- 
lation of the atmospheric response G d (z;s',s) if the atmo- 
sphere is illuminated from below with an elementary 
collimated flux beam nS= 1 in the direction s'. This dif- 
fuse Green’s function G d (z;s',s) is a solution to the radi- 
ative transfer problem adjoint to the standard radiative 
transfer equation which solution is D(z;s 0 ,s). By adjoint, 
we mean reversing the order of the atmospheric layers 
(t — > t 0 — t) but with the same boundary conditions, that 
is, G d + { 0) = 0 and G d _(z 0 ) = 0. Similar to D(z-,s 0 ,s), the 
term G d (z;s',s ) depends only on the properties of the 
atmosphere. We define the operator 

f T jL-= jdsG“{z,s,s)L\zo,s\ (A2) 

or 

which transforms Z_(t 0 ,/) into the atmospheric response. 

3. Third, the radiance field J(z;s 0 ,s) is considered as a series 
of converging terms J*\z;s 0 ,s) that quantify the radiative 
flux of photons which have undergone a number j of sur- 
face-atmosphere round trips. On the surface-atmosphere 
interface, it is possible to formalize the physical interac- 
tion that binds a term j — 1 to the next j such that 

jW(T 0 )=Rf ( t 0 ), (A3) 

and, at this level, the series converges to 

•/-(to) = I> W ( to) = [/ - R- r+] ^4%), (A4) 

J> 1 

the zeroth-order illumination of the surface being 

L+\ros) = nSe^ z °^5(s — so) + D(zo ; sos), p>0. (A5) 

[78] In all cases, the operator R expresses the reflection by 
the surface of a downward radiance field into a upward radi- 
ance field 

=R/+ ~ 1> (t 0 ) = J ds p(s' (r 0 ,s'y (A6) 

£ 2 + 

[79] Transferring the total radiance exiting the surface to 
the TOA and adding the additive contribution of the latter, 
the TOA radiation measured by the instrument becomes 

L(z = 0; So ,s) =D(0; Jo ^) + f-o, J [/-^-r+] _1 R4 1 * * * * * 0) (T 0 ), (A7) 

p < 0 . 

Al. Practical Resolution 

[so] In order to achieve the calculation, it is necessary to 
calculate each term JW(t 0 ;.s'), j> 1 using equation (A3) 
resulting in a quadruple numerical integration. In the case 
of a surface and/or atmosphere high albedo, the series can 
be long to converge and the number of terms required in 
the summation becomes important. Two main assumptions 
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become therefore necessary to circumvent the prohibitive 
computation time in this case. 

1 . After a sufficiently large number of atmosphere-surface 
round trips, two successive terms of the series become 
proportional; that is, the angular structure of the radiance 
field is preserved at the interface (t = t 0 ), the absolute 
level becoming weaker and weaker: 

J U+1) (t 0 ) = Rf+jW (t 0 ) — vJ U) (to) • (A8) 

For many natural surfaces that are reasonably anisotropic, this 
number is very small j = 2. Then we only need to calculate the 
terms j[p (t 0 ) ,j = 2,3 with equation (A3) starting with 


[82] The last step consists in transferring the radiance to 
the TOA, thanks to the direct and diffuse term of the Green’s 
function as it is done in equations ( 1 )— (4). 


Appendix B: Ross-Thick and Li-Sparse Kernels 

[83] The volumetric or Ross-Thick kernel is a single- 
scattering approximation of radiative transfer theory consist- 
ing of a Lambertian background and a layer of small scatterers 
with unifonn angle distribution and equal transmittance and 
reflectance [Roujean et al., 1992]. The form of this kernel, 
nonnalized to zero for 6q = 0°, 6 = 0°, is 


/I'^to/s) = Sg 0 e Zalllo p(so,s) 

+ i J ds'p(s ', s)p'D(t 0 ; so,s') . (A9) 

Quantity rj ~ can be now evaluated, allowing a simpli- 
fied expression of the upward radiance field requiring the 
computation of only two successive quadruple integrations 

J-{ t 0 ;s) = jW(zo;s) + J ~ . (A10) 

1 - 77 

2. Assuming that radiation is enough isotropic so that 
multiple surface-atmosphere reflections occur according 
to a Lambertian law, factor rj is the product of two bi- 
hemispherical albedos r/ = c 0 p 0 , the surface and atmo- 
spheric albedos ( p 0 and c 0 ) being 


MOoAg) 


(n/2 — g)cosg + sing n 

cosPo + cosP 4 


(Bl) 


[84] The geometric or Li-Sparse kernel assumes a sparse 
ensemble of randomly located spheroids casting shadows 
on the background, which is assumed Lambertian [ Wanner 
et al., 1995]. This geometric term is given by the proportions 
of sunlit and shaded scene components in a scene consisting 
of randomly located spheroids of height-to-center-of-crown 
h and crown vertical-to-horizontal radius ratio b/r. For 
CRISM processing, we take b/r= 1 and b/r =2 as it is done 
for MODIS processing [Lucht et al., 2000]. If the sunlit 
component is simply assumed to vary as l/cos<9 0 , the kernel 
takes on the reciprocal form 


Po=—J ds'p'J ds-p(s,s'), (All) 

a- n + 

co = -J ds'p'J dsG d (ro;s,s'). 

n+ O" 

[si] The resolution proceeds with the simplification of 
term = RT+ s) , which is the diffuse and 

direct solar illumination reflected once by the surface, then 
undergoing one surface-atmosphere round trip. As the first 
diffuse contribution is already a slowly changing function 
with respect to direction, the angular structure is conserved 
in the operation which reduces to a multiplication by 77. 
The second contribution can also be approximated but less 
drastically. This step gives the final expression for the up- 
ward radiance rising at the surface-atmosphere interface 


J-(jt Q ;s) = Sp 0 e T0 ^ 0 

1 


C0Pl(s)p 2 W 


tt[1 - C0P0] . 


p(so,s) + 

1 C0P0 
ds'p(s',s)p'D( To, ;so,s') 


(A 12) 


/g(0o, 0 , <p) = 0(00, d ’ V>) - sec 6 0 - secP 
+ - ( 1 + cosg ) sec 6 0 sec 6 , 


where 


O = - (t — sintcost) ^secP’o + secP ^ , 


cost = 


D 2 + I tanP 0 t an 0 sinp 


b secP 0 + secP 

D = \/ tan 2 P 0 + tan 2 P' — ItwiG' cos</p, 
cosg = cos P 0 cosP 4- sinP 0 sinP cos^, 

0o = tan -1 ^-tanPo^j , P = tan -1 ^-tan0^ . 


(B2) 


(B3) 


Appendix C: Rejection and Convergence Criteria 

[85] The quality of the surface solutions provided by MARS- 
ReCO is assured by the following set of rejection criteria 
adopted from Lyapustin et al. [2012]: 


where 

P\{s) = J J ds'p{s',s), p 2 (s 0 ) = J j ds'p(s 0 ,s'). (A13) 

q+ tr 


1 . Before inversion, we check if the photometric curve to be 
processed satisfies a minimum phase angle sampling, 
namely, max(g 7 ) — min(g 7 ) >20°, V/ € [ 1 ,Ng\ . Inversion 
is aborted otherwise. Additionally, MARS-ReCO is run 
only if there are at least three angular measurements as 
the inversion problem becomes unconstrained otherwise 
(note the three unknown variables k). 
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2. After retrieving the surface solution k ( " + on iteration n. 


that satisfy R c — R 


(«+l) 


> 

J 'I J J 

C ]1 pr are excluded. The inversion is restarted 


3. 


the angular measurements^ 

4 (Jrj = 4 ' 

on iteration n + 1 with a reduced number of geometries. 
Note that a R j is the a posteriori ler error bar on the mod- 
eled TO A reflectance for each geometry. Testing proved 
that a 4 a threshold is a good trade-off between rejection 
of wrong solutions and limitations of the inversion 
method. This criterion is especially useful for challeng- 
ing CRISM observations such as those acquired over 
the polar regions or under turbid atmospheric conditions. 
The physical sense of the retrieved kernel weights k <n ~ 1 * 
is verified by checking that the associated surface albedo 


q ' " 1 \d 0 ) is positive and lower than unity when 0 O = 1 5 1 


45°, 60°. The inversion of the current photometric curve 
is aborted otherwise. 


[86] At the end of each iteration, the retrieved surface so- 
lution is evaluated for convergence by the following conver- 
gence criterion: 

1. The quality of the derived surface BRF is assessed by a 
confidence index that is initially low (flag = 0, when 
n = 0). The retrieved solution k*” + 1 1 obtained on iteration 
n is compared with the previous solution k ( " ) as follows | 
q {n + l \0 o )-q w {d o )\< 0.001, where <9 0 =15°, 45°, 60°. 
Each time a new retrieval agrees with the previous solution 
according to this criterion, flag is increased by 1 and the 
inversion is repeated on next iteration n+ 1 using k (,!+1) . 
The retrieval is considered to be reliable when flag =4. This 
criterion is useful for CRISM polar observations for which 
the nonlinear term R" 1 becomes larger due to the high bright- 
ness of the surface and the extreme acquisition geometry. 

Cl. Notation 

Cl.l. Spectral and Directional Parameters 

0, 0 O : view and solar zenith angle 

p, p 0 : cosine of view and solar zenith angle 

tp, g: view-sun relative azimuth and phase angle 

s, s 0 : view and sun directions defined by (0,p) and 

(Ro.Po)- In this work, ip 0 = 0° 

A: wavelength 


C2. Atmospheric Parameters 


t 0 (2): 

aerosol optical depth 

c 0 (A): 

spherical albedo of the 
atmosphere 

G d (s 0 ,sfl): 

diffuse Green’s function of 
the atmosphere 

D(s 0 ,s,A), R D (s 0 ,s,A): 

atmospheric reflected path 
radiance and reflectance 

D s (s 0 ,s,A): 

path radiance incident on 
the surface 

F s (so, A), F d (s 0 ,A ): 

direct and diffuse radiative 
fluxes incident on the surface 

F D °""(s 0 ,A), F Up (s 0 ,A): 

downwelling and upwelling 
radiative fluxes at the surface 

L s (s 0 , s, A), Lf(s 0 ,s,F): 

direct and diffuse upwelling 
radiance after surface reflection 

a(s 0 ,i): 

multiple reflection factor 


C3. Top-of-atmosphere Parameters 

I(s 0 ,s,A)/F(s 0 ,A): ratio of measured intensity to solar flux 
Z,(.?o,s,2), R(s 0 ,s,A): gas-free radiance and reflectance 

extraterrestrial solar spectral irradiance 


C4. Surface-related Parameters 


p(s 0 ,s,A): 

q(s 0 ,A): 


f G (s 0 ,s,X),f v (s 0 ,s,A): 


k = {k\A),k G {A),k v (A)}-. 


bidirectional reflectance factor 
directional-hemispherical 
surface albedo 
RTLS geometric and 
volumetric kernels 
weights for the RTLS kernels 


C5. Inversion Parameters 


k sn i- 

{F l ,F c ',F v }\ 

R%Lo,pY 


;-(5 0 ,s,2): 


C To , C*, C r 


Fpi, C pk , Cpp. 


RMSE: 


e P : 


weights providing the best fit 
between the model and the 
CRISM measurements 
multiplicative factors for the 
kernel weights 

nonlinear term comprising the 
surface dependence of the TOA 
reflectance 

reduced measurements such that 
r = R — R d — R nl 

covariance matrix of the CRISM 
measurements, the AOD estimate, 
the RTLS weights, and the reduced 
measurements 

a posteriori covariance matrix of the 
reduced measurements, the RTLS 
weights, and the retrieved BRF 
root mean square error between the 
sensed and the modeled TOA 
reflectances 

a posteriori standard deviation of a 
retrieved photometric curve in 
BRF units 

relative error of the retrieved BRF 
photometric curve as regards the 
synthetic data 
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